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I.  INTRODUCTION 


The  calculation  of  non-adiabat ic  coupling  matrix  elements  between 
adiabatic  states  has  traditionally  been  a  very  difficult  problem.  Recent 
studies  have  relied  on  the  existence  of  a  compact  diabatic  or  pseudo-diabat ic 
basis,  or  on  finite  difference  methods.  We  have  recently  proposed  an 
analytical  method  of  obtaining  first-order  non-adiabat ic  coupling  matrix 
elements  (NACMEs) .  In  our  approach  we  employ  a  state-averaged  MCSCF  procedure 
to  define  a  common  orbital  basis.  With  a  few  simple  modifications  we  are  then 
able  to  make  use  of  the  machinery  developed  to  calculate  Cl  gradients  and 
MCSCF  second  derivatives  to  obtain  the  desired  quantities.  In  fact,  the 
methods  employed  to  obtain  MCSCF  wavefunc t ions ,  MCSCF  second  derivatives, 
multi-reference  Cl  gradients  and  NACMEs  can  be  cast  into  one  unified  and 
compact  framework. 

It  is  the  purpose  of  this  manuscript  to  outline  the  methods  used  to 
obtain  energy  derivatives  and  NACMEs.  The  first  section  will  be  devoted  to  an 
outline  of  the  second-order  MCSCF  procedure  employed  in  the  ALCHEMY  II  program 
package.  In  the  second  section,  we  show  how  MCSCF  derivatives  can  be  fit 
neatly  into  this  package  and  finally,  we  show  that  analytical  NACMEs  can  be 
obtained  within  this  framework  if  a  common  orbital  set  is  used  to  describe 
both  states. 

A.  Second-Order  MCSCF  Methodology 

In  this  section,  we  outline  the  second-order  MCSCF  method  employed  in  the 
ALCHEMY  II  series  of  programs.  The  construction  of  the  orbital  hessian  is 
completely  vectorizable  and  requires  no  formula  tape.  The  Cl  variations  are 
handled  directly  in  the  CSF  basis,  so  that  large  Cl  expansions  can  be 
addressed.  Many  of  the  matrices  constructed  as  precursors  to  the  orbital 
hessian  are  retained  and  used  to  construct  the  orbital  contribution  to  the 
inhomogenious  portion  of  the  coupled-per turbed  MCSCF  (CP-MCSCF)  equations. 
Similiarlyj  the  code  employed  to  directly  compute  the  result  of  multiplying 
the  Cl-orbi tal  coupling  portion  of  the  MCSCF  hessian  times  a  trial  vector  is 
used  to  construct  the  Cl  contribution  to  the  inhomogeneous  portion  of  the  CP- 
MCSCF  equations. 

In  the  discussions  which  follow,  I,  J,  K,  L  will  be  used  to  denote 
inactive  orbitals,  A,  B,  C,  D  will  be  used  to  denote  active  orbitals,  R,  S,  T, 
U  will  be  used  to  denote  any  orbital  in  either  the  inactive,  active,  or 
virtual  spaces,  and  a,  b,  c,  d,  will  be  used  to  denote  atomic  orbitals.  The 
wavefunction  will  be  represented  as  follows, 

♦  ■  E  Vp  (1) 

p 

where  <j>  represents  our  CSF  basis  and  Cj  is  an  element  of  the  Cl  vector.  The 
energy  expression  is  partitioned  into  an  inactive  portion,  an  inactive-active 
portion,  and  a  purely  active  portion. 


E  -  E  2hIr  *  !  (2  j”  -  K-) 
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k  and  k  are  spin-coupling  constants,  h  and  g  are  the  usual  one  and  two 

electron  Eperators  appearing  in  the  Born-Oppenheimer  hamiltonian.  D^  and 
^ABCD  are  elements  °f  the  one  and  two-particle  density  matrix 
respectively.  Unitary  orbital  variations  are  introduced  by  means  of  an 
exponential  operator  containing  only  non-redundant  terms,  ^ 


U  =  e  A  =  1-A  +  j  A2  + 
=  UQ  +  +  U2  +  . • • 


(3) 


The  elements  of  A  are  our  variational  parameters.  The 
expanded  to  second-order  in  A  and  then  differentiated, 
and  second-order  terms  arising  from  U2  are  obtained  by 
Lagrangian  matrix,^  ® 


energy  expression  is 
First-order  variations 
constructing  a 


JTA  “  DABhTB  +  2 
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CD 

^ABCD  JTB} 


(4) 


Second-order  terms  arising  from  simultaneous  first-order  variations  are 
obtained  by  first  constructing  a  Y  matrix.^- 


v  _  u  ^  ^  ,J  TCD 

Y(AT)(BS)  dAB  hTS  Z  dABCD  JTS 
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here  and  are  scaled  elements  of  the  one  and  two-particle  density 

matrix.  The  scaling  factor  depends  only  on  the  symmetry  of  the  indices  of  the 
density  matrix  elements,  (whether  A=B  ,  C=D,  or  AB=CD) .  Similiarly,  d^BCB  is 
generated  by  scaling  the  density  matrix  element.  It  is  important  to 

note  that  the  one-electron  and  coulomb  contributions  to  the  Y  matrix,  YJ,  can 
be  contracted  and  used  to  construct  the  Lagrangian.1  The  Y  matrix  and  the 

Lagrangian  are  constructed  by  multiplying  a  scaled  density  element  times  a 
coulomb  or  exchange  matrix,  a  procedure  which  is  readily  vector i zable . 


The  Y  matrix  is  partitioned  into  inactive-inactive,  inactive-active,  and 
active-active  blocks.  The  inactive-inactive  block  can  be  constructed  by 
weighting  the  coulomb  and  exchange  matrices  with  a  small  number  of 
constants.  For  example,  an  off-diagonal  matrix  in  the  inactive-inactive  sub¬ 
block  of  the  Y  matrix  is  constructed  from  only  three  terms, 


7IJ  _  TIJ  L  T  IJ  +  JI 

fTS  3IJ  JTS  +  bIJ  KTS  +  bIJ  KTS  * 


(6) 


The  inactive-active  portion  of  the  Y  matrix  requires  only  scaled  one-particle 
density  matrix  elements. 
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AI 

TS 
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(7) 


The  partitioning  of  the  Y  matrix  is  computationally  attractive  as  it  reduces 
the  number  of  matrices  which  need  to  be  held  in  core  at  one  time.  In 
addition,  it  provides  the  flexibility  of  computing  the  inactive-inactive  sub¬ 
block  or  the  inactive-active  sub-block  directly  as  discussed  by  Bacskay^  and 
Olsen,  et  al . ,  respectively.  In  addition,  the  occupied-virtual,  occupied- 
virtual  portions  of  the  matrix  can  be  combined  with  Lagrangian  contributions 
to  generate  an  occupied-virtual,  occupied-virtual  sub-blocks  of  the  orbital 
hessian,  which  need  not  be  kept  in  core. 


The  construction  of  the  orbital  hessian  proceeds  as  follows, 


(1)  A  reduced  integral  transformation  is  used  to  form  the  J  and  K 
matrices  in  the  ao  basis  and  to  calculate  the  integrals  needed  to  perform  the 
Cl, 


(2)  In  a  two  step  MCSCF  procedure  the  Cl  calculation  is  then  carried 

out , 


(3)  The  Cl  vector  is  used  to  construct  the  unique  elements  of  the  one 
and  two-particle  density  matrix  over  the  active  orbitals, 


(4)  The  density  matrix  elements  are  scaled  and  reordered  in  the  case  of 
exchange  contributions, 
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(5)  The  coulomb  and  fock  operator  contributions  to  the  Y  matrix  are 
generated , 

(6)  The  Y  matrix  is  contracted  with  the  molecular  orbitals  to  produce 
the  active  portion  of  the  Lagrangian.  The  inactive  portion  of  the  Lagrangian 
is  generated  from  fock  matrices, 

(7)  The  exchange  contributions  are  added  to  the  Y  matrix,  and 

(8)  The  Lagrangian  is  folded  to  produce  the  MCSCF  orbital  gradient  (g ^ 

=  -  Lg^)  anc*  the  Y  matrix  is  transformed  to  the  mo  basis,  similarly  folded 

and  combined  with  Lagrangian  contributions  to  produce  the  MCSCF  hessian. 


The  Newt on-Raphson  linear  equations  are  solved  iteratively, 


where  60  and  6C  represent  orbital  and  Cl  variations,  respectively,  and  the 
remaining  portions  of  the  MCSCF  hessian  are  treated  directly.^  ^  The 
second-order  Cl  terms  are  handled  with  a  direct-CI  program. 

Multiplication  of  the  Cl-orbital  portion  of  the  hessian  times  a  trial 
vector  are  re-expressed  as  a  gradient  constructed  with  transition  density 
matrices,  1 


r  9  E  ^  op  _  <$C 

303 c  6C  80 


(9) 


or  updated  integrals, 


12 


4o>  «0  -  -  (H-E)S0|c>  ■ 


(10) 


The  updated  integrals  used  to  construct  g  are  obtained  as  follows, 

<ABCD>  60  =  E  (60  <TBCD>  +  60  <ATCD>  (11) 

_  TA  TB 

T 

+  60m„  <ABTD>  +  60__  <ABCT>) 

TC  TD 

where  <TBCD>  =  <TB|g|CD>.  The  integrals  with  three  indices  transformed 
<TBCD>,  etc.,  have  been  generated  earlier  by  contracting  a  coulomb  matrix. 
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II.  MCSCF  SECOND-DERIVATIVES 


The  formalism  needed  to  generate  a  MCSCF  force  constant  matrix  has  been 
discussed  by  a  number  of  authors,  j  13-16  ^he  presented  by  Page,  et 

al . ,  has  been  incorporated  in  the  MCSCF  framework  discussed  above  as  it 
seems  to  offer  a  number  of  computational  advantages.  In  particular,  the 
construction  of  the  inhoraogeneity  in  the  CP-MCSCF  equations, 


(12) 


involves  the  same  code  employed  in  the  construction  of  the  MCSCF  Lagrangian 
and  in  the  handling  of  the  Cl-orbital  coupling  terms  in  the  direct  solution  to 
the  Newton-Raphson  equations.  In  this  formalism,  ga  and  g01  are  the  orbital 
ajjtd  Cl  gradients  constructed  with  derivative  ao  integrals.  Similiarly, 
gQ  and  gc  are  gradients  constructed  with  updated  integrals  where  one  index 
has  been  transformed  with  the  Ta  matrix.  Ta  is  an  upper  triangular  matrix 
constructed  from  transforming  the  derivative  ao  overlap  integrals  into  the  mo 
basis , 


S>R 


2 


(13) 


and 

Tsr  =0,  S>R 
Ta 

g  is  generated  by  contracting  the  Y  matrix  with  Ta  and  transforming  the 
MCSCF  Lagrangian  with  Ta ,  to  first  produce  the  Ta  Lagrangian,  LTa  . 


I 

MS 


AM 

RS 


+  T.  T 
S 


a 

SR 


JSA 


(14) 


where  the  summation  over  M  runs  over  all  occupied  orbitals.  The  code  employed 
in  the  MCSCF  program  to  produce  the  Lagrangian  from  is  also  used  at  this 
point  to  contract  the  Y  matrix. 

Ta  . 

gc  is  obtained  from  the  code  emplgged  in  the  direct  solution  of  the 
Newton-Raphson  equations  to  construct  gp  .  In  this  case  the  updated  integrals 
are  obtained  by  contracting  the  partially  transformed  integrals,  <TBCD>,  etc., 
with  the  Ta  matrix. 


The  final  expression  for  the  first  and  second  derivatives  of  the  MCSCF 
energy  expressions  are,  ^ 
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where  Up  is  the  derivative  of  the  mo  coefficient  matrix  obtained  by  combining 
the  solution  of  the  CP-MCSCF  equations  and  T^  .  These  equations  only  employ 
quantities  (La  and  L  a)  used  to  construct  the  inhomogeneity  in  the  CP-MCSCF 
equations,  the  solution  to  the  CP-MCSCF  equations  themselves  or  terms 
generated  from  a  SCF  derivative  package  after  the  MCSCF  density  matrix  has 
been  transformed  to  the  ao  basis.  The  actual  assemblage  of  the  MCSCF  force- 
constant  matrix,  after  the  solution  the  CP-MCSCF  equations  and  the  contraction 
of  the  derivative  ao  integrals  with  ao  density  matrix,  is  extremely  simple. 


III.  NON- ADIABATIC  COUPLING  MATRIX  ELEMENTS 
The  first-order  non-adiabat ic  coupling  matrix  element, (NACME) , 


d,m  ..  \  = 

(M,N,a)  M'da  N  M 1  N 


(17) 


can  be  calculated  in  a  straightforward  manner  if  the  same  orbital  basis  is 
used  to  represent  the  two  states. ^  ^  In  this  case,  the  required 
wavefunction  derivatives  are  obtained  from  a  state-averaged  CP-MCSCF  equation. 
The  NACME  is  broken  down  into  two  terms, 


D  /  N  —  D  +  D 

(M,N3ct)  C  0 


=  <cM  ca>  +  z  dm’n  (ua  +  va  )* 

Ml  N  ‘  ST  ST  ST 
ST 


(18) 


The  first  term  is  simply  the  overlap  of  a  Cl  vector  and  a  derivative  Cl 

vector.  For  a  MCSCF  wavefunction,  C?T  is  obtained  from  the  solution  to  the  CP- 

N 

MCSCF  equations,  so  this  term  poses  no  new  problems.  For  a  Cl  wavefunction, 
the  direct  evaluation  of  this  term  would  require  the  solution  of  at  least  one 
CP-CI  type  of  equation, 
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(19) 


<»-%>  c“  - 


-  <h-V(“,cm 


Multiplying  Eq.  19  on  the  left  by  CN,  we  find  the  usual  perturbation 
expression,  ^ 


<CJC“>  =  (E  -E„) 


-1 


jni^m 


M  N 


<cnIh 


(a) 


V 


(20) 


The  term  on  the  right-hand-side  of  this  equation  is  a  familiar  one.  It  is 
similar  to  the  expression  for  a  Cl  gradient.  The  only  difference  is  that  we 
need  to  employ  transition  density  matrices  in  place  of  normal  density  matrices 
in  our  equations.  The  formal  expression  for  a  Cl  gradient  is** 


F(«) 

EM 
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ab  ab 
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abed 


D 


M,M 

abed 


Dabcd 


Z 

TS 


rM,M 


TS  TS 


(21) 
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The  first  term  is  the  trace  of  the  density  matrix  with  derivative  ao 
integrals,  while  the  second  term  is  the  contraction  of  the  Cl  Lagrangian  with 
the  derivative  of  the  mo  coefficients.  Thus,^ 


Dc  = 


nM’N  uCt  A  ^M,N  a 

,  ab  ab  abed  gabcd 

ab  abed 


+  E  L^’N  l£  }  .  (22) 

TS  lb  Ib 


Therefore ,  we  able  to  make  use  of  the  same  code  used  to  compute  Cl  gradients. 

In  the  expression  for  the  orbital  contribution  to  the  first-order  NACME 
we  are  dealing  with  a  one-electron  operator  so  the  resulting  expression 
involves  the  trace  of  an  one-particle  density  matrix  and  an  overlap  term,^ 


Do  = 


ST 


dm’n  (ua  +  va  ) 

ST  s  T  ST7 


(23) 


Dg’T  is  an  element  of  a  non- symmetric ,  one-particle  transition  density 
mafrix.  Ug  T  is  the  derivative  of  the  mo  coefficients  obtained  from  the 
solution  to* the  state-average  CP-MCSCF  equations,  and  T  is  obtained  by 
transforming  half-derivative  ao  overlap  integrals  into  the  mo  basis.  Two 
overlap  terms  are  needed  as  the  derivative  of  the  molecular  orbitals  are 
expressed  as,^® 


dO 

da 


’  3^  (xt>  ’ 


x  t  +  xt 
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=  xat  +  OU 
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(24) 
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where  0  are  the  molecular  orbitals,  x  is  a  vector  of  atomic  orbitals,  and  t  is 
the  mo  coefficient  matrix.  ^ 


ab 

and  <xa|x^>  is  a  half-derivative  overlap  integral. 

The  new  terms  required  to  compute  a  NACME  are  a  square  one-particle 
transition  matrix  and  half-derivative  ao  overlap  integrals.  The  Cl  gradient 
package  must  also  be  modified  to  produce  the  required  one  and  two-particle 
transition  density  matrices,  while  the  MCSCF  package  need  only  be  modified  to 
compute  a  square  one-particle  transition  density  matrix. 


c  <x 

aS  a 


(25) 


This  formalism  can  be  extended  to  analytically  compute  second-order 
NACMEs.  This  requires  the  set-up  of  the  second-order  state-averaged  CP-MCSCF 
equations  and  the  code  to  compute  a  square  two-particle  transition  density 
matrix.  A  second-order  Cl  NACME  would  require  the  solution  of  the  first-order 
CP-CI  equations. 


IV.  CONCLUSION 

A  unified  treatment  of  energy  derivatives  and  non-adiabat ic  coupling 
matrix  elements  has  been  outlined.  We  showed  that  in  the  case  of  both  energy 
derivatives  and  NACMEs  the  final  equations  could  be  rewritten  in  a  form  that 
resulted  in  the  efficient  use  of  existing  MCSCF  or  derivative  methods. 
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